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Major  Contributions  Include: 

1)  All-atom  computational  chemistry  (molecular  dynamics  -  MD)  simulation  of  normal  shock  waves  is 
now  possible.  Pure  MD  simulations  were  validated  with  experimental  data  for  shock  structure.  The  only 
model  input  into  such  simulations  is  the  potential  energy  surface  (PES)  that  governs  individual  atomic 
interaction  forces,  developed  by  chemists  and  physicists. 

2)  In  contrast  to  existing  rotational  energy  models,  we  found  that  the  rotational  relaxation  rate  (Trot)  depends 
strongly  on  the  initial  degree  of  nonequilibrium  and  the  direction  towards  the  equilibrium  state. 
Compressing  flows  involve  fast  excitation  whereas  expanding  flows  involve  slow  relaxation  for  the  same 
equilibrium  temperature.  We  developed  a  new  rotation  model  for  both  DSMC  and  CFD  simulation. 

3)  At  high  temperatures  we  found  clear  evidence  of  rotational- vibrational  energy  coupling.  Existing  models 
using  relaxation  times  (Trot  and  Tvib)  are  not  able  to  reproduce  the  ro-vibrational  coupling  found  in  MD 
simulations.  We  propose  that  an  additional  relaxation  time  (Tr<,t-vib)  is  required  to  model  direct  energy 
transfer  between  rotation  and  vibration  at  high  temperatures  and  we  presented  a  new  ro-vibrational 
model.  Preliminary  simulations  for  dissociating  nitrogen  were  also  demonstrated. 

4)  A  numerical  method  that  replaces  the  collision  model  in  DSMC  with  MD  trajectories  was  developed 
and  shown  to  reproduce  exactly  pure  MD  results.  This  enables  simulation  of  axisymmetric  and  even  3D 
flows  where  a  potential  energy  surface  is  the  only  model  input.  Thus,  the  accuracy  of  pure  MD  is 
maintained  for  full  flow  fields;  directly  linking  computational  chemistry  with  aerothermodynamics. 
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II.  Introduction  and  Existing  State-of-the-Art 

Currently,  the  two  main  computational  tools  used  by  the  aerothermodynamics  community  to  model 
hypersonic  flows  are:  (1)  Computational  Fluid  Dynamics  (CFD)  methods,  which  solve  the  multi¬ 
temperature  Navier-Stokes  equations  with  chemical  reaction  source  terms  and,  (2)  the  direct  simulation 
Monte  Carlo  (DSMC)  particle  method  [1],  which  simulates  the  Boltzmann  equation  [2]  and  is  therefore 
accurate  for  highly  nonequilibrium  flows  relevant  to  rarefied  flows  and  sharp  flow  features  with  small 
length  scales.  Currently,  both  CFD  and  DSMC  use  essentially  the  same  physical  models  for  rotational- 
vibrational  excitation  and  dissociation  phenomenon,  which  are  based  on  a  limited  number  of  near¬ 
equilibrium  experiments  performed  at  low  temperatures.  The  purpose  of  the  current  research  is  to  develop 
a  new  modeling  capability,  based  on  computational  chemistry,  to  provide  a  more  fundamental 
understanding  and  develop  more  accurate  thermochemical  models  for  CFD  and  DSMC. 

Despite  operating  on  molecular  quantities,  DSMC  currently  employs  the  same  empirical  transport, 
internal  energy,  and  chemical  models  used  in  CFD  methods.  For  example,  the  models  most  widely 
employed  in  CFD  codes  include  the  Park  model  for  air  dissociation  [3],  the  Millikan  and  White  rates  for 
vibrational  relaxation  [4],  Blottner  curve  fits  for  viscosity  [5]  and  Wilke’s  mixing  rule  for  diffusion  [6]. 
These  rates  are  parameterized  as  being  temperature  dependent  (typically  in  Arhenius  or  power-law  form). 
For  DSMC,  these  empirically  constructed  transport,  relaxation,  and  reaction  rates  are  simply  converted 
from  functions  of  temperature  and  re-written  as  probabilities  based  on  relative  collision  energy,  then 
subsequently  applied  to  individual  DSMC  particle  collisions  [7-9].  Thus  under  continuum  gas  conditions, 
DSMC  results  are  ensured  to  be  consistent  with  CFD  results  (by  their  very  construction).  However,  this 
Top-down’  (phenomenological)  modeling  approach  is  inherently  limiting,  since  empirical  rates  for 
complex  interactions  are  often  absent  (difficult  to  determine  experimentally)  or  greatly  extrapolated  for 
hypersonic  conditions.  Furthermore,  their  accuracy  in  non-equilibrium  regions  is  certainly  questionable. 

For  example,  discrepancies  between  high-speed  wind-tunnel  experimental  measurements  and  numerical 
simulations  are  attributed  to  the  uncertainty  in  vibrational  relaxation  and  dissociation  rates  in  high 
temperature  air  [10].  Very  recently,  discrepancy  between  such  experiments  and  simulations  has  been 
attributed  to  possible  order-of-magnitude  uncertainty  in  oxygen  recombination  rates  [11].  Furthermore, 
the  lack  of  collision  cross-section  data  leads  to  significant  uncertainty  in  diffusive  transport  through  high- 
temperature  boundary  layers  [12],  and  therefore  uncertainties  in  predicted  surface  heating  rates.  Finally, 
collisional  non-equilibrium  has  been  shown  to  be  significant  even  at  altitudes  where  continuum  (CFD) 
analysis  is  thought  to  be  accurate  [7]  and  may  therefore  require  non-equilibrium  internal  energy/chemistry 
modeling  and  kinetic  simulation  (DSMC). 

The  uncertainty  in  such  rates  is  a  major  hurdle  facing  predictive  simulation  capabilities  for  hypersonic 
flows.  This  has  prompted  the  development  of  extended  vibrational  models  [13]  and  even  models  where 
each  rotational  and  vibrational  energy  level  is  represented  as  its  own  species  in  computationally  intensive 
“state-to-state”  simulations.  Recently,  computational  chemistry  researchers  at  NASA  Ames  Research 
Center  have  used  quantum  computations  to  fit  very  detailed  Potential  Energy  Surfaces  (PES)  for 
interactions  between  two  and  three  nitrogen  atoms,  and  are  currently  working  on  PES  for  four  atoms 
[14,15].  These  quantum-computed  PES,  corresponding  to  all  possible  inter-atomic  orientations,  are  used 
to  ‘fit’  a  force-field  that  can  be  used  in  so-called  Quasi-Classical-Trajectory  (QCT)  Molecular  Dynamics 
(MD)  simulations  of  collisions.  NASA  researchers  have  simulated  millions  of  individual  collisions  using 
QCT-MD  and  compiled  a  database  of  collision  cross-sections  for  the  purpose  of  integration  to  obtain 
temperature-dependent  rate  data  for  use  in  CFD  models  [15].  For  use  in  the  DSMC  method,  Motsumoto 
et  al  developed  the  Dynamic  Molecular  Collision  (DMC)  model  [16]  specifically  for  the  translational- 
rotational  energy  transfer  in  diatomic  nitrogen.  Many  isolated  collisions  were  performed  using  MD  with 
an  appropriate  interatomic  potential  and  a  database  of  collision  cross-sections  and  pre/post  energy 
distributions  was  compiled  and  used  within  the  DSMC  method.  More  recently,  instead  of  performing 
isolated  MD  collisions,  Kotov  and  Surzhikov  [17]  and  Kantor  et  al.  [18]  simulated  the  dissociation  of  air, 
and  oxygen-hydrogen  combustion,  respectively,  in  a  periodic  box  using  pure  MD  for  the  entire  system  of 
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molecules.  However,  both  studies  report  being  severely  limited  in  the  number  of  simulated  molecules, 
being  restricted  to  very  high  densities.  The  proposed  research  seeks  to  vastly  improve  upon  such  research. 

The  goal  of  the  current  research  is  to  develop  complete  4bottom-up’  gas-phase  collision  models  for  high 
temperature  dissociated  air  for  use  in  both  DSMC  and  CFD  methods.  A  key  advancement  required  to 
achieve  this  goal  is  the  development  of  accelerated  computational  chemistry  (molecular  dynamics  -MD) 
methods  designed  specifically  for  dilute  gases.  The  accelerated  MD  methods  developed  by  our  research 
are  capable  of  directly  incorporating  state-of-the-art  interatomic  potentials  (that  govern  the  interaction 
forces  between  individual  atoms)  from  computational  chemistry  as  the  only  modeling  assumption.  Such  a 
capability  enables  new  DSMC  and  CFD  model  development  and  a  deeper  understanding  of  energy 
exchange  and  reaction  processes  in  hypersonic  flows  at  the  most  fundamental  level.  In  the  original 
proposal  we  described  ideas  to  develop  pure  MD  simulations  and  accelerated  MD  simulations  that  enable 
the  simulation  of  ID  hypersonic  flow  features  (such  as  shock  waves  and  stagnation  line  profiles). 
However,  our  research  has  progressed  dramatically  beyond  this,  and  we  are  able  to  demonstrate 
accelerated  MD  simulations  of  2D  and  even  3D  flow  fields  where  only  the  interactions  between  atoms  are 
considered  (no  collision  or  rate  models  are  required).  The  result  of  our  research  is  that  future  advances 
from  computational  chemistry  can  be  directly  incorporated  into  large-scale  simulations  of  full  3D 
hypersonic  flow  fields  in  strong  thermochemical  nonequilibrium. 
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III.  The  Molecular  Gas  Dynamic  Simulator  (MGDS)  DSMC  Code 

My  research  group’s  core  capability  is  large-scale  stochastic  molecular  simulation;  specifically  the  direct 
simulation  Monte  Carlo  (DSMC)  method.  With  partial  support  from  this  grant,  we  have  developed  a  new 
state-of-the-art,  parallel  DSMC  code  called  the  Molecular  Gas  Dynamic  Simulator  (MGDS)  code.  The 
code  has  a  number  of  novel  features.  MGDS  uses  an  embedded  three-level  Cartesian  flow  grid  with 
automated  adaptive  mesh  refinement  (AMR).  The  refinement  is  arbitrary  (non-binary)  and  enables 
accurate  and  efficient  simulations  with  little  user  input  [19-21].  In  addition,  MGDS  contains  a  robust 
“cut-cell”  subroutine  that  cuts  complex  3D  geometry  from  the  background  Cartesian  grid  and  exactly 
computes  the  volumes  of  all  cut  cells  [22].  Combined 
within  a  DSMC  solver,  these  two  features  enable 
molecular-level  physics  to  be  applied  to  real 
engineering  problems.  For  example,  Fig.  1  shows 
hypersonic  flow  over  complex  satellite  geometry  with 
sharp  leading  edges  and  Fig.  2  shows  the  solution  and 
adapted  flow  field  grid  for  hypersonic  flow  over  a 
blunt  body.  It  is  noted  that  since  DSMC  is  a  particle 
method,  the  flow  field  grid  does  not  need  to  be  aligned 
with  sharp  gradients  as  required  by  CFD  methods, 
which  solve  stiff  sets  of  discretized  partial  differential 
equations.  This  enables  accurate  DSMC  simulations  to 
be  performed  for  flows  with  sharp  gradients  on 
complex  3D  geometry  with  very  little  user- time 
required  (no  gridding). 

In  addition  to  the  practicality  for  complex  3D  flows,  Figure  1:  DSMC  simulation  of  a  satellite  geometry, 
our  new  DSMC  code  acts  as  a  bridge  between 

computational  chemistry  modeling  and  continuum  fluid  mechanics.  With  existing  parallel  computer 
clusters,  DSMC  is  able  to  perform  near-continuum  simulations  using  only  molecular  physics  models.  I 
believe  such  research  has  the  potential  to  make  a  large  impact  on  the  field  of  CFD.  Upcoming  sections 
will  describe  pure  computational  chemistry  simulations  and  section  VII  will  describe  how  the 
research  performed  under  this  grant  has  led  to  a  new  method;  one  which  embeds  computational 
chemistry  directly  into  large-scale  DSMC  simulations,  so  that  simulations  such  as  those  in  Figs.  1 
and  2  can  now  be  performed  using  only  atomic  interaction  forces  as  the  model  input. 


Figure  2:  DSMC  solution  for  a  blunt-body  flow  requiring  significant  mesh  adaptation  in  the  fore-body  and  wake 
regions  (left).  Close  up  image  of  the  solution  and  3 -level  Cartesian  flow  field  grid  near  the  shoulder  region  (right). 
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IV.  All-Atom  Molecular  Dynamics  of  Shock  Waves  in  Noble  Gases 

Through  this  grant,  my  research  group  has  demonstrated  that  pure  computational  chemistry  of  simulation 
of  shock  waves  and  shock  layers  is  computationally  feasible  with  modest  parallel  computing  resources. 
Such  simulations,  which  are  herein  referred  to  as  “all-atom  molecular  dynamics”  (MD)  simulations, 
require  only  a  potential  energy  surface  (called  a  “PES”)  that  dictates  the  interaction  forces  between 
individual  atoms,  as  the  model  input.  Thus,  unlike  DSMC  or  CFD,  no  models  for  viscosity,  diffusion, 
thermal  conductivity,  internal  energy  relaxation,  chemical  reactions,  cross-sections,  or  even  state-resolved 
probabilities/rates  are  required.  Rather,  every  real  atom  in  the  system  (such  as  the  thin  column  containing 
a  normal  shock  depicted  in  Fig.  3)  is  simulated.  Each  atom  interacts  with  its  neighboring  atoms  through  a 
potential  energy  surface  (PES),  such  as  the  simple  Lennard-Jones  (LJ)  PES  [23]  shown  in  the  inset  of  Fig. 
3,  which  is  given  by  the  following  equation: 


ipinj)  =  4e 
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where  r/;-  is  the  distance  between  atoms  i  and  /,  and  8  and  a  are  constants  specific  to  the  species  (/  and  j)  of 
the  interacting  atoms.  The  motion  of  all  atoms  is  then  computed  by  integrating  Newton’s  law  of  motion: 

n=vh 


N 


F  i{r a)  =  mvi  =  2 


j=  i 


dipin,)  ru 

dra  ru 


where  N  is  the  total  number  of  neighboring  atoms  that  influence  atom  i  (usually  prescribed  to  be  within  a 
conservative  cutoff  radius  of  atom  i).  Thus,  such  simulations  are  fundamental  and  extremely  general,  as 
the  only  simulation  input  is  the  PES  (for  example,  the  LJ  PES  shown  above).  The  highly  scalable  MD 
code,  LAMMPS  [24],  can  be  used  to  perform  the  simultaneous  integration  of  many  millions  of  atoms. 
LAMMPS  is  freely  available  from  Sandia  National  Laboratory  [25]  and  thus  such  MD  simulations  are 
repeatable  by  other  researchers  in  the  field. 
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IX 


Figure  3:  Schematic  of  an  all-atom  molecular  dynamics  (MD)  simulation  of  a  normal  shock  wave.  Inset  shows  a 
simple  Lennard-Jones  (LJ)  potential  energy  surface  (PES)  dictating  atomic  interaction  forces. 


The  main  point  of  this  section  is  to  demonstrate  that  for  gas  systems  where  accurate  PES  have  been  well- 
established,  such  pure  MD  simulations  of  shock  waves  can  be  compared  directly  with  experimental 
results.  We  have  carefully  verified  and  validated  our  pure  MD  simulations  with  a  large  amount  of 
experimental  data.  Figure  4  shows  results  for  pure  argon,  where  the  shock  thickness  computed  by  MD  is 
in  excellent  agreement  with  experimental  data  from  Alsmeyer  and  the  simulations  are  able  to  precisely 
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resolve  the  bi-modal  velocity  distribution  functions  throughout  the  shock  wave.  Complete  details  are 
fully  described  in  our  recent  journal  articles  [26,27]. 


Figure  4:  Pure  MD  simulation  results  for  normal  shock  waves  in  argon.  Shock  thickness  results  vs.  experiment 
(left)  and  velocity  distribution  functions  within  the  shock  (right). 


We  then  performed  pure  MD  simulations  of  mixtures  of  noble  gases  for  which  a  great  deal  of 
experimental  data  exists.  Results  for  argon-helium  mixtures  are  shown  in  Figs.  5  and  6,  and  results  of 
xenon-helium  mixtures  are  shown  in  Fig.  7.  For  such  flows,  there  is  a  large  mass  disparity  between 
helium  atoms  (light)  compared  to  argon  and  xenon  atoms  (heavy).  This  causes  a  species  concentration 
separation  within  the  shock  wave  due  to  diffusion  effects,  and  the  simulations  are  in  excellent  agreement 
with  experimental  results.  These  pure  MD  simulations  also  demonstrate  that  very  precise  results  are 
computationally  feasible  even  for  mixtures  with  trace  species  (as  low  as  1.5%  in  some  simulations). 
Complete  details  are  fully  described  in  our  recent  journal  article  [28]. 


Figure  5:  Argon  temperature  profiles  within  a  normal  shock  wave  in  an  argon-helium  mixture.  Symbols  are 
experimental  data  (7*:  triangles,  Ty=Tz:  circles),  lines  are  MD  simulation  results. 
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Figure  6:  Species  velocity  profiles  within  a  normal  shock  wave  in  an  argon-helium  mixture.  Symbols  are 
experimental  data  (Argon  :  triangles,  Helium :  circles),  lines  are  MD  simulation  results. 


Figure  7:  Species  temperature  profiles  within  a  normal  shock  wave  in  an  xenon-helium  mixture.  Symbols  are 
experimental  data  (Helium:  triangles,  Xenon :  circles),  lines  are  MD  simulation  results. 


V.  All- Atom  Molecular  Dynamics  of  Shock  Waves  in  Nitrogen  (Rotational  Relaxation) 

After  carefully  validating  our  MD  simulations  for  monatomic  gas  mixtures,  we  moved  on  to  diatomic 
gases  where  rotational  excitation  of  the  molecules  was  now  an  important  physical  aspect.  We  performed 
simulations  of  shock  waves  in  nitrogen,  using  a  site-to-site  LJ  potential,  that  is  accurate  for  nitrogen  at 
relatively  low  temperatures  (<3000  K)  where  vibrational  excitation  can  be  safely  ignored.  We  first 
validate  the  PES  by  computing  the  viscosity  of  nitrogen  over  a  wide  temperature  range  and 
comparing  with  experimental  data.  Validating  the  viscosity  (essentially  the  collision  rate)  predicted 
by  a  PES  is  a  crucial  step  that  should  always  be  performed  first.  More  advanced  physics,  such  as 
rotational  and  vibrational  relaxation,  and  chemical  reactions  are  always  defined  relative  to  the  collision 
rate.  For  example,  “collision  numbers”  (Zrot,  Zvib)  and  reaction  probabilities  are  used  as  models  that 
determine  what  fraction  of  collisions  results  in  internal  energy  exchange  and  reactions.  Thus,  if  the 
collision  rate  (i.e.  the  viscosity)  is  not  modeled  correctly,  then  the  interpreted  rates  of  internal  energy 
transfer  and  chemical  reactions  will  be  incorrect  as  well.  We  determine  the  viscosity  predicted  by  a  PES 
using  the  Green-Kubo  theory  [29,30],  i.e.,  by  computing  the  autocorrelation  function  of  the  off-diagonal 
pressure  tensor  terms,  which  were  estimated  using  the  virial  equation.  After  validating  the  viscosity 
prediction,  we  performed  normal  shock  calculations  and  compared  the  results  to  experimental  data  from 
Robben  and  Talbot.  Figure  8  shows  MD  results  compared  to  experiment  for  the  density  and  rotational 
temperature  profiles  in  Mach  7  and  12.7  nitrogen  shock  waves.  The  experimental  data  set  also  included 
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rotational  energy  distribution  functions  throughout  the  shock.  These  are  plotted  in  in  Fig.  8,  compared  to 
MD  simulation  results.  The  distributions  are  seen  to  start  at  a  Boltzmann  Distribution  (BD)  corresponding 
to  the  upstream  temperature  and  relaxing  to  a  BD  corresponding  to  the  downstream  temperature. 
However,  within  the  shock  wave,  the  rotational  energy  distributions  are  seen  to  be  non-Boltzmann. 

Complete  details  are  fully  described  in  our  recent  journal  article  [31]. 


P  (Exp.) 
Tro,  (Exp.) 
P  (MD) 


Figure  8:  MD  simulation  results  compared  to  experiment  for  nitrogen  shock  waves.  Density  and  rotational 
temperature  profiles  within  the  shock  waves  (left).  Rotational  energy  distribution  functions  (right). 
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After  validating  our  nitrogen  simulations  with 
experimental  data,  we  were  able  to  analyze  the  large 
amount  of  data  produced  by  such  pure  MD  simulations  ^ 

to  solve  a  long-standing  problem  regarding  P 

translational-rotational  energy  transfer.  Specifically,  it  ^ 
had  been  postulated  in  the  literature  that  the  rate  of  P 
rotational-translation  energy  transfer  in  a  gas  may  not  be  a  ^ 
function  of  the  equilibrium  temperature  only  (the  approach 
taken  by  all  existing  models),  but  instead,  may  depend  T 
strongly  on  the  initial  degree  of  nonequilibrium  and  the 
direction  towards  the  equilibrium  state  (i.e.  compressing 
vs.  expanding  flows).  In  a  1998  review  article  on  this  topic 
[32],  Wysong  and  Wadsworth  concluded  that,  “ the 
primary  complicating  factor ,  and  that  which  may  be  of 
most  use  in  improving  phenomenological  models  of  the 
process ,  is  the  potentially  large  effect  of  the  size  and  sign  of 
the  initial  deviation  from  equilibrium ”.  In  addition  to  normal 
shocks,  we  performed  isothermal  relaxation  calculations  with 

pure  MD  where  molecules  were  initialized  (in  a  periodic  box  -  zero-dimensional)  with  their  rotational 
energies  out  of  equilibrium  with  translation  and  they  relaxed  via  collisions  to  the  equilibrium  state.  As 
seen  in  Fig.  9,  we  observed  that,  for  the  same  equilibrium  temperature,  the  rate  of  rotational  energy 
relaxation  varied  significantly  depending  on  the  initial  state.  Specifically,  our  MD  simulations  clearly 
showed  that  rotational  relaxation  rates  were  a  strong  function  of  both  the  direction  to  equilibrium  and  the 
initial  degree  (magnitude)  of  nonequilibrium.  Figure  10  shows  how,  for  the  same  equilibrium  temperature 
(Tinf),  compressing  flows  (Ttr0  -  Trot°>  0)  have  short  relaxation  times  (Trot),  while  for  expanding  flows  (Ttr0 
-  Trot0  <  0)  relaxation  times  are  much  larger. 


2E-09 


Figure  9:  Rotaional  temperature  profiles  for 
isothermal  relaxation  MD  simulations  for 
different  initial  conditions. 


Currently,  the  most  widely  used  rotation  model  is  the  Parker  model,  which  only  models  rotational 
relaxation  (xrot  or  equivalently  Zrot)  in  function  of  the  equilibrium  temperature.  Therefore,  the  Parker 
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model  only  accounts  for  a  line  (Ttr0-Trot0  =  0)  extracted  from  our 
data  in  Fig.  10.  The  Parker  model  along  with  our  data  extracted  on 
this  line  are  shown  in  Fig.  11  along  with  other  equilibrium  results 
from  computational  chemistry.  Also  plotted  in  Fig.  11  is  the 
available  experimental  data  that  is  seen  to  have  a  huge  variation, 
and  was  only  interpreted/parameterized  in  terms  of  the  equilibrium 
temperature.  Our  results,  shown  fully  in  Fig.  10,  are  much  more 
general  and  conclude  that  the  magnitude  of  initial  nonequilibrium 
and  direction  to  the  equilibrium  state  have  a  much  stronger  effect 
on  rotational  relaxation  that  the  equilibrium  temperature.  Enabled 
by  our  MD  simulation  capability,  we  were  the  first  to  quantify 
these  effects  for  nitrogen.  Our  most  recent  paper,  in  the  AIAA 
Journal  [33],  presents  a  new  “direction-dependent”  rotation  model 
(the  NDD  model)  for  use  in  both  DSMC  and  CFD.  In  addition  to 
the  new  model,  this  article  makes  a  number  of  further  contributions 
including  a  mathematical  formulation  that  ensures  detailed  balance 
for  DSMC  models,  and  a  mathematical  framework  that  rigorously 
integrates  a  DSMC  model  (in  the  limit  of  equilibrium)  to  obtain  a 
consistent  model  for  CFD.  The  new  DSMC  model  parameterizes 
the  probability  of  rotational  energy  exchange  as  a  function  of  both 
the  rotational  energy  and  translational  energy  involved  in  a 
collision: 


«—  Tjnf=  100  K  (ISOTHERMAL) 

• -  Tjnf  =  300  K  (ISOTHERMAL) 

Tjnf=  1000  K  (ISOTHERMAL) 


Tjnf  =  2000  K  (ISOTHERMAL) 


Figure  10:  Rotational  relaxation  rates 
predicted  by  MD  simulation  in  function 
of  both  Ttr  and  Trot. 
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When  integrated  in  the  limit  of  equilibrium  (Maxwell-Boltzmann  translational  and  rotational  energy 
distribution  functions),  the  following  temperature  dependent  model  results,  which  is  suitable  for  use  in 
multi-temperature  Navier-Stokes  (CFD)  simulations: 


ZrotiTuTr)  — 


Finally,  this  new  NDD  model  was 
used  to  compute  full  shock  wave 
profiles  and  compare  to  pure  MD 
solutions.  We  found  that  the  new 
NDD  model,  which  correctly  models 
the  direction-dependence  of  the 
relaxation  rate,  is  able  to  accurately 
reproduce  full  MD  results  while  the 
existing  Parker  model  is  not  (Fig.  12). 
Furthermore,  the  NDD  model  was 
also  able  to  successfully  predict 
rotational  relaxation  in  expanding 
flows,  while  the  Parker  model  was 
not.  Thus  this  research,  used  novel 
pure  MD  simulations  to  discover  new 


Molecular  Dynamics  (this  work) 

Parker 

Semi-classical  (Billing  and  Wang,  1992) 
Classical  rigid  rotor  (Nyeland  and  Billing,  1988) 
Annis  and  Malinauskas  (Exp.,  1971) 

Carnevale  et  al.  (Exp.,  1967) 

Ganzi  and  Sandler  (Exp.,  1971) 

Healy  and  Stonick  (Exp.,  1969) 
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Figure  11:  Rotational  relaxation  rate  in  the  limit  of  thermal  equilibrium 
(Ttr=Trot)  from  various  sources. 
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physical  processes  import  in  rotational  relaxation.  These  high  fidelity  calculations  provided  the  data 
needed  to  create  new  DSMC  and  CFD  models.  Complete  details  are  fully  described  in  our  recent 
journal  articles  [31,33]. 


x  [gm] 

(a)  NDD  model 

Figure  12:  Normal  shock  calculations  using  DSMC  wi 
compared  with  a  solution  from  pure  MD  simulation. 
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(b)  Parker  model 

the  new  NDD  model  (a)  and  the  Parker  model  (b) 


VI.  All-Atom  Molecular  Dynamics  of  Shock  Waves  in  Nitrogen  (Rotational-Vibrational  Excitation) 

After  studying  translational-rotational  energy  exchange,  we  moved  on  to  study  translational-rotational- 
vibrational  energy  exchange  and  dissociation.  The  temperature  range  considered  is  now  much  higher  and 
it  is  known  that  vibrational  processes  are  sensitive  to  the  “repulsive  wall”  region  of  the  PES  (very  close 
range  atomic  interactions).  In  this  energy  range,  we  use  a  potential  energy  surface  (PES)  from  Ling  and 
Rigby  (LR)  [34]  that  has  been  used  by  a  number  of  researchers  to  study  rotation  and  vibration  in  nitrogen 
at  high  temperatures  [35,36].  The  LR  PES  that  governs  the  interaction  between  an  atom  and  atoms  of 
neighboring  molecules  (the  mtermolecular  potential)  is  given  by: 

c 

0(ry)  =  Dexp(-arij  -  /9r?  )  -  /d(ry  )-^ 

rij 


«nJ)  =  exP(-i[(^)2-i]) 


and  is  parameterized  for  diatomic  nitrogen  [34].  The  interactions  between  atoms  of  the  same  molecule 
(the  m^ramolecular  potential)  is  described  by  the  Harmonic  Oscillator  potential  (HO),  parameterized  for 
nitrogen.  This  means  that  for  simulations  using  the  LR  +  HO  potentials  that  molecules  cannot  dissociate 
and  only  rotational-vibrational  processes  are  simulated.  In  later  simulations  we  use  increasing  degrees  of 
anharmonicity  and  finally  a  Morse-type  potential  (anharmonic  oscillator  -  AHO),  which  allows  for 
molecular  bond  breaking  and  dissociation.  It  is  important  to  note  that  simulations  performed  with  any 
of  these  PES  operate  on  individual  atoms  only  (position  and  velocity  of  atoms).  Thus  there  is  no 
assumed  decoupling  of  rotational  and  vibrational  energies  within  a  simulation.  Only  as  a  post¬ 
processing  step  do  we  analyze  the  positions  and  velocities  of  atoms  bonded  in  nitrogen  molecules  to 
compute  rotation  and  vibrational  energies.  This  is  a  unique  aspect  of  such  MD  simulations  and  enables 
us  to  study  ro-vibrational  coupling  effects  where  no  other  numerical  method  can. 
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We  performed  a  large  number  of  isothermal  relaxation  calculations,  similar  to  those  performed  for  Fig.  9, 
except  now  at  much  higher  temperatures,  where  significant  vibrational  energy  exchange  occurs.  Two 
sample  results  are  shown  in  Fig.  13,  and  the  results  are  very  interesting. 


t  [ns] 


Figure  13:  Isothermal  relaxation  simulations  using  pure  MD  with  rotating  and  vibrating  nitrogen.  Symbols  are 
MD  solutions  and  lines  are  solutions  to  the  modified  Jeans  and  Landau-Teller  equations  that  we  propose. 

In  Fig.  13  (a),  initially  Trot=Tvib  «  Tt  and  as  the  MD  simulation  proceeds  the  rotational  energy  increases 
first,  followed  by  the  vibrational  temperature.  In  Fig.  13  (b),  initially  Tvib  «  Trot  =  Tt  and  as  the 
simulation  proceeds  the  vibrational  energy  increases,  however,  the  rotational  energy  actually  decreases 
first,  before  increasing  back  into  equilibrium  with  translation.  Both  of  these  initial  conditions  are 
representative  of  the  gas  state  immediately  downstream  of  a  strong  shock  wave.  We  find  that  no  existing 
rotational  or  vibrational  energy  transfer  models  are  capable  of  matching  these  trends. 

At  the  continuum  level,  the  equation  typically  used  to  describe  the  relaxation  process  is 

of  the  form: 

de^  _ 

dt  n  U 

where  s*(Tt )  =  (i/2  k #Ti(£ )  represents  the  instantaneous  equilibrium  energy  associated  with 

the  i-th  internal  energy  mode  with  0  degrees  of  freedom  and  ks  is  Boltzmann’s  constant. 

Tt(t)  is  the  time-dependent  translational  temperature  of  the  system.  For  rotation  (vibration), 

Eq.  2  is  referred  to  as  the  Jeans  (Landau- Teller)  equation.  It  is  common  to  express 

Ti  =  ZiTc,  (2) 

where  r;  is  the  mean  collision  time  (a  function  of  Tt )  and  Z.t  is  a  so-called  relaxation  collision 
number. 

For  an  isothermal  relaxation,  where  the  translational  temperature  (Tt)  remains  a  constant,  the  solution  to 
equation  (1)  for  the  rotational  temperature  profile  is: 

4C(I)  =  AT®,  exp(— =  A  7*  exp  {-tvc/Zrol) 
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For  lower  temperatures  we  found  that  all  MD  rotational  relaxation  solutions  were  well-fit  by  such  an 
expression  (for  example  see  the  curves  in  Fig.  9  and  the  Zrot  values  in  Fig.  10).  However,  the  solution  for 
Trot  in  Fig.  13  (a)  cannot  be  fit  by  such  an  expression.  Likewise,  since  equation  (1)  relaxes  either 
rotational  (or  vibrational)  energy  into  equilibrium  with  only  translational  energy  (which  is  fixed  for  the 
simulations  in  Fig..  13),  the  decrease  in  rotational  temperature  seen  in  Fig.  13  (b)  cannot  be  predicted  by 
equation  (1).  What  is  happening  at  these  high  temperatures  is  that  there  is  energy  exchange  directly 
between  rotational  and  vibrational  energy  modes  during  collisions  of  nitrogen  molecules.  We  find  that  an 
additional  term  is  required  in  Eq.  (1)  in  order  to  account  for  this  physical  process.  Specifically: 

lCT,(t)  -  »( rTr(t)  t  0  _  |C -  lC,T,(t)  n 

7r  +u  dt  -  Tv 

Q  _  |(Cr  +  ~  |(Cr  +  C v)Tr{t) 

Trv 


That  is,  in  order  to  match  the  relaxation  behavior  seen  in  our  MD  calculations  one  must  consider 
not  only  xv  and  xr  (relaxation  defined  in  reference  to  translational  energy),  but  also  a  coupled 
relaxation  parameter,  xrv.  We  were  then  able  to  determine  each  of  the  parameters  xv  ,  xr ,  and  x^  for 
nitrogen  from  our  MD  simulations;  the  results  are  shown  in  Fig.  14.  In  order  to  interpret  all  parameters 
from  our  MD  simulations  we  had  to  first  determine  the  viscosity  (^xc)  of  the  LR  +  HO  PES  over  the  entire 
temperature  range,  followed  by  determining  the  direction-dependent  rotational  relaxation  behavior  of  this 
potential  (xrot  -  very  similar  to  that  presented  in  Fig.  10,  but  over  a  much  larger  temperature  range), 
followed  by  a  validation  of  the  predicted  vibrational  relaxation  at  low  temperatures  compared  to  the 
experimental  data  of  Millikan  and  White.  This,  then  allowed  us  to  determine  the  x^  required  to  account 
for  the  rotational  energy  coupling  to  vibration  (as  seen  in  Fig.  13  (b)),  and  finally  enabled  us  to  determine 
the  resulting  xvib  to  match  all  MD  profiles. 
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Ling-Rigby+HO  (MD,  present  work) 
Ling-Rigby+HO  (CTC  DSMC,  present  work) 
Millikan-White  (MW)  correlation  (Exp.,  1963) 
MW,  high-temperature  correction  (Park,  1985) 
MW,  high-temperature  correction  (Boyd,  1999) 
SCTC  (Billing  and  Fisher,  1979) 

SCTC  (POT  I,  Cacciatore  et  al.,  2005) 
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Figure  14:  Relaxation  time  constants  and  collision  numbers  (xr  -  Zr,  Tv  -  ZV|  and  Trv  —  Zrv)  determined  by  pure  MD 
simulations  in  function  of  equilibrium  (translational)  temperature. 
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We  report  two  major  findings.  First,  rotational  and  vibrational  relaxation  times  (Zvib  and  Zrot) 
approach  each  other  at  high  temperatures  (Fig.  14  -  right),  and  thus  both  processes  are  equally 
important.  Whereas,  existing  CFD  models  account  only  for  vibrational  relaxation  and  assume  that 
rotational  energy  is  in  equilibrium  with  translational  energy.  Second,  in  order  to  match  the  relaxation 
behavior  seen  in  our  MD  calculations  one  must  consider  not  only  xv  and  %r  (relaxation  defined  in 
reference  to  translational  energy),  but  also  a  coupled  relaxation  parameter,  xrv.  The  magnitude  of  the 
new  term  (which  includes  x^),  is  such  that  it  only  becomes  significant  at  high  temperatures,  while  at  low 
temperatures,  the  original  Jeans  and  Landau-Teller  equations  are  recovered.  The  additional  term 
represents  relaxation  of  rotational  and  vibrational  energy  modes  to  an  equilibrium  internal  energy  in 
addition  to  the  original  equations,  which  represent  relaxation  to  an  equilibrium  translational  energy.  Our 
proposed  modification  alters  the  energy  transfer  between  energy  modes  at  high  temperatures  where  these 
processes  are  coupled  to  dissociation.  Complete  details  are  fully  described  in  a  journal  article  we  have 
just  submitted  for  publication  and  that  is  currently  under  review  [37]. 

Although  the  temperatures  in  Figs.  13  and  14 
are  well  above  the  equilibrium  dissociation 
temperature  for  N2,  we  note  that  in  a 
nonequilibrium  shock  layer  there  will  be 
substantial  molecular  collisions  within  this 
range  of  translational  energies  before  the  gas- 
molecules  are  fully  dissociated.  Furthermore, 
although  not  shown  in  this  report,  our 
preliminary  results  show  that  increasing  the 
anharmonicity  of  the  PES  (more  realistic  of 
actual  N2  molecules  that  dissociate)  shifts  the 
rotational-vibration  coupling  effects  to  lower 
temperatures.  Figure  15  shows  a  preliminary 
pure  MD  simulation  of  a  nitrogen  shock 
wave  including  dissociation.  This  simulation 
includes  trans-rot-vib  energy  transfer  with 
ro-vibrational  coupling  and  dissociation. 

We  note  that  the  scatter  is  due  to  the  absence 
of  atomic  species  in  the  upstream  portion  and 
the  absence  of  diatomics  the  downstream 
portion. 

It  is  not  yet  clear  to  what  degree  such 
rotational-vibrational  coupling  ultimately  affects  dissociation.  However,  the  unique  modeling  capabilities 
and  expertise  we  have  developed  due  to  this  AFOSR  grant  enable  my  research  group  to  unravel  these 
physical  processes  and  understand  high  speed  flows  at  the  most  fundamental  level.  We  are  working  to 
confirm  and  refine  these  findings  using  new  state-of-the-art  PES  that  have  recently  been  constructed  as 
part  of  a  large  Multi-University  Research  Initiative  (MURI)  centered  in  the  Aerospace  Engineering  & 
Mechanics  and  Chemistry  departments  at  Minnesota.  If  we  can  confirm  our  findings  and  quantify  the 
affect  on  dissociation  chemistry,  this  would  be  a  very  significant  contribution  to  the  high-speed  flow 
community. 


Figure  15:  Preliminary  pure  MD  simulation  of  a  nitrogen  shock 
wave  including  dissociation. 


VII.  Embedding  MD  trajectories  within  DSMC  Simulations 

Finally,  it  is  possible  to  imbed  molecular  dynamics  trajectories  directly  into  DSMC  simulations.  In  this 
manner,  the  rigorous  simplifications  employed  by  DSMC  are  maintained;  using  simulator  particles  that 
each  represent  a  large  number  of  identical  real  molecules,  moving  simulator  particles  with  timesteps  on 
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the  order  of  xc,  and  stochastically  selecting  collision  pairs  and  initial  orientations  within  volumes 
(computational  cells)  on  the  order  of  X.  However,  the  collision  models  in  DSMC  are  completely  replaced 
by  trajectory  calculations  performed  on  arbitrary  PES.  This  technique  was  introduced  by  Koura  in  2001  as 
Classical-Trajectory-Calculation  (CTC)  DSMC  [38-40]  where  trajectories  were  combined  with  the  null- 
collision  DSMC  algorithm.  As  part  of  the  current  grant,  we  published  a  CTC-DSMC  implementation  that 
uses  the  No-Time-Counter  (NTC)  scheme  (the  most  popular  DSMC  scheme)  with  a  cross-section  that  is 
determined  by  the  PES  and  is  consistently  applied  as  the  maximum  impact  parameter  for  the  trajectory 
calculations  [41].  Since  both  DSMC-CTC  and  pure  MD  simulations  can  use  the  same  PES  as  the  only 
model  input,  the  two  particle  simulation  techniques  can  now  be  directly  compared.  For  example,  Fig.  16 
compares  CTC-DSMC  and  pure  MD  solutions  for  a  Mach  7  diatomic  nitrogen  shock  where  a  comparison 
between  normalized  density  and  temperature  profiles  is  shown.  The  results  of  CTC-DSMC  and  MD 
simulations  are  in  perfect  agreement.  The  rotational  energy  distribution  functions  at  different  values  of  the 
normalized  Trot  within  the  shock  are  also  found  to  be  exactly  equivalent.  The  MD  density  and  rotational 
temperature  profiles  and  rotational  distribution  functions  have  been  previously  validated  with 
experimental  data  [31].  It  is  noted  that  no  DSMC  collision  model  exists  that  is  able  to  exactly  reproduce 
the  MD  results  to  this  level  of  precision  [31,33].  Thus  for  the  flows  considered  (and  potentially  for  more 
complex  flows),  CTC-DSMC  is  shown  to  be  purely  a  numerical  acceleration  technique  for  the  MD 
simulation  of  dilute  gases,  since  both  methods  require  only  the  PES  as  a  model  input.  Complete  details 
are  fully  described  in  our  recent  journal  article  [41]. 

Although  much  more  efficient  than  pure  MD  simulation  (simulating  every  real  atom  with  femtosecond 
timescales),  CTC-DSMC  remains  roughly  2  to  3  orders  of  magnitude  more  computationally  expensive 
than  standard  DSMC,  depending  on  the  complexity  of  the  PES  used  for  the  trajectories.  However,  we 
have  shown  that  CTC-DSMC  calculations  scale  extremely  well  on  parallel  architectures  (including  GPU 
architectures  [41]).  Essentially,  the  CTC-DSMC  technique  automates  the  generation  of  state-to-state 
transition  rate  data.  Instead  of 
performing  billions  of  independent 
trajectories  to  determine  state- 
transition  probabilities,  the  billions  of 
trajectories  are  imbedded  within  an 
actual  flow  simulation.  In  this  manner, 
the  trajectories  performed  within  a 
CTC-DSMC  simulation  are  naturally 
those  that  are  most  relevant  to  the 
problem  being  simulated.  The 
influence  of  specific  PES  can 
immediately  be  determined  and 
variations  in  PES  (which  model  the 
same  species  interactions  to  different 
accuracy)  can  be  directly  compared 
for  full  flow  fields  without  the  need  to 
generate  and  converge  full  state-to- 
state  cross-section  models  and 
implement  in  a  DSMC  code  for  each 
variation  of  a  PES.  Thus,  CTC-DSMC 
may  be  useful  in  guiding  the 
development  of  state-to-state  models 
and  lead  to  insight  into  physical-based 
model  reduction  strategies. 


O  MD 


Figure  16:  Comparison  between  CTC-DSMC  and  pure  MD  solutions  for  a 
shock  wave  in  diatomic  nitrogen. 
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Finally,  we  have  demonstrated  that  CTC- 
DSMC  capability  allows  solutions  to  full 
3D  nonequilibrium  flow  fields  over 
complex  geometry,  with  accuracy  that  is 
expected  to  be  equivalent  to  full  MD 
simulation.  An  example  3D  solution  is 
shown  in  Fig.  17  below  for  hypersonic 
flow  of  argon  over  a  capsule  geometry. 
Although  this  is  only  argon  gas,  the 
important  fact  is  that  this  shock-layer 
solution  was  obtained  with  only  a  single 
PES  as  input  to  the  simulation,  and  no 
other  model.  Future  work  will  continue  to 
develop  this  capability  and  incorporate 
state-of-the-art  PES.  While  such 
simulations  are  not  expected  to  overlap 
with  the  3D  near-continuum  flows  in  the 
near  future,  they  certainly  overlap  with 
rarefied  DSMC  simulations  and  for 
2D/axi- symmetric  flows,  may  overlap 
with  CFD  solutions  if  thousands  of 
processors  are  used  (which  is  feasible 
even  now). 


Figure  17:  Three-dimensional  DSMC-CTC  solution  for  hypersonic 
flow  of  argon  over  a  capsule  geometry.  This  simulation  required 
only  128  cores  for  12  hours  using  the  MGDS  code. 


VIII.  Conclusions 

1 .  All-atom  computational  chemistry  (molecular  dynamics  -  MD)  simulation  of  normal  shock  waves  is 
possible.  A  simulation  methodology  was  developed  and  rigorously  validated  with  experimental  data  for 
shock  structure  in  mixtures  of  nobles  gases  and  diatomic  nitrogen.  Precise  MD  simulations,  able  to 
resolve  distribution  functions  and  also  species  with  low  concentrations,  are  possible  with  moderate 
computational  resources  (100  core  processors  for  a  few  days). 

2.  All-atom  MD  simulations  of  nitrogen  discovered  clear  “direction-dependence”  of  translational- 
rotational  energy  transfer  that  is  not  captured  by  existing  models.  We  found  that  compressing  flows 
involve  fast  rotational  excitation,  whereas  expanding  flows  involve  slower  relaxation,  for  the  same 
equilibrium  temperature.  We  developed  a  new  model  for  both  DSMC  and  CFD  that  reproduces 
experimental  data  and  is  consistent  between  MD,  DSMC,  and  CFD. 

3.  All-atom  MD  simulations  of  nitrogen  discovered  rotation- vibration  coupling  at  high  temperatures.  It 
was  found  that  the  widely  used  Jeans  and  Landau-Teller  equations  for  rotation  and  vibrational  energy 
transfer  (in  the  multi-temperature  Navier-Stokes  equations)  are  unable  to  reproduce  the  ro-vibrational 
coupling  found  in  MD  simulations.  We  propose  that  an  additional  relaxation  time  (Trot-vib)  is  required  to 
model  direct  energy  transfer  between  rotation  and  vibration  at  high  temperatures  and  present  a  new  ro- 
vibrational  model.  Preliminary  simulations  of  dissociating  nitrogen  are  also  demonstrated  and  future  work 
will  use  new,  state-of-the-art,  PES  being  developed  for  air  species  interactions  at  NASA  Ames  Research 
Center  and  also  through  a  large  AFOSR  MURI  grant  centered  at  the  University  of  Minnesota. 

4.  We  developed  a  combined  MD-DSMC  technique  that  replaces  the  collision  model  in  DSMC  with  MD 
trajectories  and  verified  the  method  to  reproduce  exactly  pure  MD  results.  This  is  a  significant 
advancement  that  enables  simulation  of  axisymmetric  and  even  3D  flows  where  a  potential  energy  surface 
is  the  only  model  input.  Thus,  the  accuracy  of  pure  MD  is  maintained  for  full  flow  fields;  directly  linking 
computational  chemistry  with  aerothermodynamics. 
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5.  In  general,  this  grant  has  developed  expertise  in  computational  chemistry  techniques  themselves,  their 
validation,  their  application  to  high  temperature  flows  of  interest,  and  how  to  use  their  predictions  to  form 
reduced  order  models  for  DSMC  and  CFD.  Future  work  can  combine  all  of  these  developments  with  new 
state-of-the-art  PES  being  developed  by  computational  chemists  with  the  goal  of  developing  new  models 
that  are  consistent  across  all  scales  (MD  -  DSMC  -  CFD  -  Experiment).  Such  research  would  enable  an 
understanding  of  hypersonic  flows  at  the  most  fundamental  level. 
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